library(car)
library(lsmeans)
library(calibrate)
library(dplyr)
#library(DEGreport)
library(DESeq2)
library(DEFormats)
library(edgeR)
library(emmeans)
library(ggpubr)
library(ggsci)
library(ggplot2)
library(gridExtra)
library(pheatmap)
library(reshape2)
library(RColorBrewer)
library(scales)
library(rstatix)
library(tidyr)
library(magrittr)
library(PCAtools)
library(tidyverse)
library(Biobase)
library(marray)
library(limma)
library(gplots)
data <- read.csv("../../SerumAntibodies_2021/IgG_MCF_SH_179_raw.csv", header=T, row.names = 1)
#Columns 16-47 are empty and are therefore removed
data <- as.data.frame((data))[,-c(16:47)]
#Creating metadata including Strain and age
Strain <- c(rep(c("BALBc", "NOD"), each=6), rep(c("NOD"), each=3))
Age <- c(rep(c("8wks", "28wks"), each=3), rep(c("8wks", "24wks", "33wks"), each=3))
Group <- c(rep(c("A","B","C","D","E"), each=3))
colData <- as.data.frame(cbind(c('B1_8', 'B2_8', 'B3_8', 'B1_28', 'B2_28', 'B3_28', 'N1_8', 'N2_8', 'N3_8', 'N1_24', 'N2_24', 'N3_24', 'N1_33', 'N2_33', 'N3_33'), Strain, Age, Group))
colnames(colData) <- c('Sample', "Strain", "Age", "Group")
rownames(colData) <- colData$Sample
colnames(data) <- colData$Sample
colData$Strain <- factor(colData$Strain)
colData$Strain <- relevel(colData$Strain, ref = "BALBc")
colData$Age <- factor(colData$Age, levels = c("8wks", "24wks", "28wks", "33wks"))
colData$Age <- relevel(colData$Age, ref='8wks')
#Boxplot of raw Signal Distribution, not normalized
mydata<-as.matrix(log2(data+0.5))
boxplot(as.data.frame(mydata),main="Signal distribution, Not normalized",xlab="Samples",ylab="Log2 Intensity",cex=0.8)
######2 Clustering samples based on Pearson correlation ############
matcor<-cor(mydata)
range(as.vector(matcor))
## [1] 0.8503868 1.0000000
strain<-as.character(as.numeric(colData$Strain))
age<-as.character(as.numeric(colData$Age))
library(gplots)
heatmap.2(matcor,trace="none",col=heat.colors(40),ColSideColors=strain,
RowSideColors=age, cexCol=1,cexRow=1,labCol="")
mydata <- as.matrix((data))
dataN <- (mydata[1:95,]/mydata[96,] * 100) #df with all rows
dataN <- mydata
boxplot(as.data.frame(log2(dataN+0.5)),main="Signal distribution, Not normalized",xlab="Samples",ylab="Log2 Intensity",cex=0.8)
library(matrixStats) #for function rowmedians
IgG_SNR <- read.csv("../../SerumAntibodies_2021/IgG_MCF_SNR.csv", header=T, row.names = 1)[,1:15]
IgG_SNR$average <- rowMeans(IgG_SNR)
IgG_SNR$med <- rowMedians(as.matrix(IgG_SNR))
mydata <- data[which(IgG_SNR$average>3),] #filtering rows with high noise (i.e., rows with Signal to Noise Ratio or SNR > 3)
mydata<-as.matrix(mydata)
dataN <- mydata
boxplot(as.data.frame(log2(dataN+0.5)),main="Signal distribution, Not normalized SNR>3",xlab="Samples",ylab="Log2 Intensity",cex=0.8)
dataN <- (mydata[1:69,]/mydata[70,]) * 100000
boxplot(as.data.frame(log10(dataN+0.5)),main="IgG normalization",col=strain)
library(limma) #for function plotMDS
plotMDS(log2(dataN+0.5))
mm <- model.matrix(~0 + Group)
y <- voom((dataN), mm, plot = T)
countData = as.data.frame(dataN)
df_dseq = melt(countData, variable.name = "Samples", value.name = "count") # reshape the matrix
df_dseq = data.frame(df_dseq, Condition = substr(df_dseq$Samples,1,1))
mycolors <- colorRampPalette(brewer.pal(8, "Set1"))(13)
ggplot(df_dseq, aes(x = log10(count+0.5), colour = Samples, fill = Samples)) +
geom_density(alpha = 0, size = 0.8) +
facet_wrap(~Condition, ncol=3) #+
#theme_minimal() + xlim(-1.5,6) +
#scale_colour_manual(values=mycolors, name="") + guides(fill="none")
matcor<-cor(dataN)
range(as.vector(matcor))
## [1] 0.1740427 1.0000000
strain<-as.character(as.numeric(colData$Strain))
age<-as.character(as.numeric(colData$Age))
heatmap.2(matcor,trace="none",col=heat.colors(40),ColSideColors=strain,RowSideColors=age, cexCol=1,cexRow=1,labCol="")
#mydata <- as.matrix((data+0.5))
conditions<- paste(colData$Strain[1:12], colData$Age[1:12],sep=".")
conditions <- factor(conditions, levels=unique(conditions))
design <- model.matrix(~0+ conditions)
colnames(design) <- levels(conditions)
#contrasts of the previous analysis.
#cont.matrix<- makeContrasts(B8vsN8 = NOD.8wks - BALBc.8wks, B28vN24 = NOD.24wks - BALBc.28wks,
# B28vN33 = NOD.33wks - BALBc.28wks,
# N24vN8 = NOD.24wks - NOD.8wks,
# N33vN8 = NOD.33wks - NOD.8wks,
# N33vN24 = NOD.33wks - NOD.24wks,
# levels = design)
#Introducing new contrasts by removing the 33 week old samples.
fit <- lmFit(dataN[,1:12], design)
cont.matrix<- makeContrasts(
B8vsN8 = NOD.8wks - BALBc.8wks,
B28vN24 = NOD.24wks - BALBc.28wks,
N24vN8 = NOD.24wks - NOD.8wks,
B28vB8 = BALBc.28wks - BALBc.8wks,
levels = design)
fit.cont<- contrasts.fit(fit, cont.matrix)
fit.cont<- eBayes(fit.cont)
decide <- matrix(c("fdr",0.05,
"fdr",0.1, "none",0.005, "none", 0.01),nrow=4,ncol=2,byr=T)
mysum <- as.list(1:nrow(decide))
mynum <- 0
maxmax <- 0
for (test in 1:nrow(decide)){
results<-decideTests(fit.cont,
adjust.method=decide[test,1],
p=as.numeric(decide[test,2]))
summary(results) -> mysum[[test]]
mynum[test] <-length(which(apply(results,1,function(x)any(x,na.rm=T))))
maxmax <- max(c(maxmax, as.vector(mysum[[test]][c(1,3),])))
}
par(mfrow=c(1,nrow(decide)))
for (test in 1:nrow(decide))
{
as.numeric(as.vector(mysum[[test]][3,]))->plotMe1
as.numeric(as.vector(mysum[[test]][1,]))->plotMe2
maxData = max(plotMe1)
maxData2 = max(plotMe2)
barplot(plotMe1,horiz=T,col="red",xlim=c(-maxmax,maxmax),
main=paste("Gene Changes \np<",decide[test,2], ", " , decide[test,1],
" (" ,mynum[test] ,")",sep=""))->yy
barplot(-plotMe2,horiz=T,col="green",add=T)->yy
xx<-vector("integer",ncol(mysum[[test]]))
text(xx,yy,colnames(mysum[[test]]))
text((plotMe1+10)*0 + .9*maxData,yy+0.1,format(plotMe1,digits=3))
text((-plotMe2-10)*0 - .9*maxData2,yy+0.1,format(plotMe2,digits=3))
}
results<-decideTests(fit.cont,adjust.method="none", p=0.01)
summary(results)
## B8vsN8 B28vN24 N24vN8 B28vB8
## Down 0 0 0 0
## NotSig 68 65 67 69
## Up 1 4 2 0
write.fit(fit.cont,file="IgG_DGE_0423.csv",adjust="none",results=results, sep=',')
fitObj <- read.csv("IgG_DGE_0423.csv", header=T, row.names = 1)
toptable1 <- topTable(fit.cont, coef=2, adjust.method="none", p=0.01, sort.by = "P", n = Inf)
write.csv(toptable1,file="IgG_DGE_0423_II.csv")
myNames<-names(fitObj)
res.col<- which(regexpr("Res.",myNames)>0)
DElist<- which(apply(fitObj[,res.col],1,function(x)any(x,na.rm=T)))
length(DElist)
## [1] 5
fitObjDE <-fitObj[DElist,]
fitObjDE$Autoantigen <- rownames(fitObjDE)
fitObjDE
## AveExpr Coef.B8vsN8 Coef.B28vN24 Coef.N24vN8 Coef.B28vB8
## Collagen VI 344353.360 306988.513 698156.44 466614.955 75447.0273
## Myosin 34663.933 68677.811 42702.87 -21566.851 4408.0887
## PCNA 25915.298 16328.886 79655.24 64814.316 1487.9658
## Sm 6144.959 6763.382 14148.64 7357.128 -28.1265
## ssDNA 753270.453 679724.733 1795054.96 1226257.802 110927.5744
## t.B8vsN8 t.B28vN24 t.N24vN8 t.B28vB8 P.value.B8vsN8
## Collagen VI 2.6252712 5.970419 3.990347 0.645199725 0.028906317
## Myosin 4.6927709 2.917897 -1.473668 0.301205734 0.001322884
## PCNA 0.9255585 4.515041 3.673823 0.084341295 0.380216625
## Sm 1.7581652 3.677988 1.912512 -0.007311583 0.114578304
## ssDNA 1.4040943 3.708010 2.533057 0.229140955 0.195776818
## P.value.B28vN24 P.value.N24vN8 P.value.B28vB8 F
## Collagen VI 0.000264043 0.003544981 0.5358308 17.760699
## Myosin 0.018153085 0.176632442 0.7704991 10.407852
## PCNA 0.001686093 0.005656377 0.9347390 9.434717
## Sm 0.005621207 0.090039827 0.9943351 6.144546
## ssDNA 0.005374459 0.033497880 0.8241801 6.511897
## F.p.value Results.B8vsN8 Results.B28vN24 Results.N24vN8
## Collagen VI 0.0005205418 0 1 1
## Myosin 0.0032778110 1 0 0
## PCNA 0.0044953318 0 1 1
## Sm 0.0162032926 0 1 0
## ssDNA 0.0137576330 0 1 0
## Results.B28vB8 Autoantigen
## Collagen VI 0 Collagen VI
## Myosin 0 Myosin
## PCNA 0 PCNA
## Sm 0 Sm
## ssDNA 0 ssDNA
library(devtools)
#install_github("dpgaile/AutoAntArrayExmpl")
library(AutoAntArrayExmpl)
library(devtools)
#install_github("dpgaile/AutoAntArrayExmpl")
library(AutoAntArrayExmpl)
library(NMF)
library(quantreg)
library(asbio)
library(fdrtool)
#library(discreteMTP)
#Sample Spot and Background Distributions
IgG_NSI <- read.csv("../../SerumAntibodies_2021/IgG_MCF_NSI.csv", header=T, row.names = 1)[1:96,1:12]
IgG_SNR <- read.csv("../../SerumAntibodies_2021/IgG_MCF_SNR.csv", header=T, row.names = 1)[1:96,1:12]
Strain <- c(rep(c("BALBc", "NOD"), each=6))
Age <- c(rep(c("8wks", "28wks"), each=3), rep(c("8wks", "24wks"), each=3))
Group <- c(rep(c("A","B","C","D"), each=3))
colData <- as.data.frame(cbind(c('B1_8', 'B2_8', 'B3_8', 'B1_28', 'B2_28', 'B3_28', 'N1_8', 'N2_8', 'N3_8', 'N1_24', 'N2_24', 'N3_24'), Strain, Age, Group))
colnames(colData) <- c('Sample', "Strain", "Age", "Group")
rownames(colData) <- colData$Sample
colnames(IgG_NSI) <- colData$Sample
colnames(IgG_SNR) <- colData$Sample
colData$Strain <- factor(colData$Strain)
colData$Strain <- relevel(colData$Strain, ref = "BALBc")
colData$Age <- factor(colData$Age)
colData$Age <- relevel(colData$Age, ref='8wks')
######Data Filtering #######
colData <- colData[1:12,] #remove the last group - M NOD 33 wks
IgG_SNR$average <- rowMeans(IgG_SNR)
IgG_SNR$med <- rowMedians(as.matrix(IgG_SNR))
IgG_NSI <- IgG_NSI[which(IgG_SNR$average>3),] #filtering rows with high noise (i.e., rows with SNR > 3)
IgG_SNR <- IgG_SNR[which(IgG_SNR$average>3),]
IgG_raw=list()
IgG_raw$NSI <- as.matrix(IgG_NSI)
IgG_raw$SNR <- as.matrix(IgG_SNR)[,1:12]
IgG_raw$SInfo <- colData
clrs=c(rep(pal_jco("default")(5), each=3))[1:12]
pchs=c(rep(3,3),rep(4,3), rep(5,3), rep(6,3), rep(7,3))[1:12]
#tiff("fig1_.tiff", units="in", width=7, height=7.4, res=300)
par(mfrow=c(2,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
plot(as.vector(IgG_raw$NSI),as.vector(IgG_raw$SNR),type="n",ylab="IgG Background Intensity (SNR)",xlab="IgG Spot Intensity (NSI)" )
abline(0,1,lwd=3,col="steelblue")
for(j in 1:12) points(IgG_raw$NSI[,j],col=clrs[j],IgG_raw$SNR[,j], pch=j,cex=0.5)
plot(log(as.vector(IgG_raw$NSI)),
log(as.vector(IgG_raw$SNR)),
type="n",ylab="log IgG Background Intensity (SNR)",xlab="log IgG Spot Intensity (NSI)")
abline(0,1, lwd=3,col="steelblue")
for(j in 1:12) points(log(IgG_raw$NSI[,j]),log(IgG_raw$SNR[,j]),
pch=j,col=clrs[j],cex=0.5)
legend(6.5,0,lty=1,lwd=3,col=clrs[c(1,4,7,11)],pch=pchs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)
plot(density(log10(IgG_raw$NSI[,1]+ 0.5),type="l"), main=" log10(NSI +0.5)", sub=" . ")
for(j in 2:12) {
dens<-density(log10(IgG_raw$NSI[,j]+0.5))
lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
legend(-0.8,0.6,lty=1,lwd=3,col=clrs[c(1,4,7,11)],pch=pchs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)
plot(density(log10(IgG_raw$NSI[,j]+0.5)), main="log10(SNR +0.5) ", sub=" . ", type="l")
for(j in 1:12) {
dens<-density(log10(IgG_raw$NSI[,j]+0.5))
lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
#dev.off()
#Centering and Scaling Differences Across Features
# Tukey Tri-mean for location
TriMnG=RowTriMeans(IgG_raw$NSI)
# Biweight midvariance for spread
bw.estG=unlist(apply(IgG_raw$NSI,1,r.bw))
#tiff("fig2.tiff", units="in", width=7, height=3.6, res=300)
par(mfrow=c(1,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,3.25,1.5,0.75))
plot(TriMnG,bw.estG,pch=c(rep(12)),
col=c(rep("darkgreen",12)),
xlab="Autoantigen Tukey TriMean",
ylab="Autoantigen Biweight Midvariance")
legend("topleft",pch=c(12),col=c("darkgreen"),
legend="IgG",bty="n")
plot(log(TriMnG),log(bw.estG),pch=c(rep(12)),
col=c(rep("midnightblue",12)),
xlab="Log Autoantigen Tukey TriMean",
ylab="Log Autoantigen Biweight Midvariance")
legend("topleft",pch=c(12),col=c("midnightblue"),
legend="IgG",bty="n")
#dev.off()
#Heatmaps heatmaps using the uniform rank-its (i.e., rank values scaled to the unit interval) of the assays values across samples (i.e., the rank-its of the spot quantifications for each autoantigen across all 10 samples)
if(!file.exists("docs/Fheatmap3.png")){
nmf.options(grid.patch=TRUE)
RnkSmplW=matrix(nrow=70,ncol=12)
RnkSmplW[,1]=IgG_raw$NSI[,1]
for(j in 2:12){ RnkSmplW[,j]=IgG_raw$NSI[,j]}
for(i in 1:70) {RnkSmplW[i,]=rank(RnkSmplW[i,])/(13)} #(row ranks)
colnames(RnkSmplW)=as.character((colData$Sample))
rownames(RnkSmplW) <- rownames(IgG_raw$NSI)
sink("mess.txt")
aheatmap(RnkSmplW,annCol=factor(colData$Strain),
main="Rank-it (Across Samples) Raw Signal IgG",
dist="euclidean",hclustfun="ward.D",
filename="Fheatmap3.png")
sink()
}
library(pheatmap)
#tiff("fig3_.tiff", units="in", width=6, height=12.5, res=300)
pheatmap(RnkSmplW, annCol=factor(colData$Strain), main="Rank-it (Across Samples) Raw Signal IgG", dist="euclidean", hclust="ward")
#dev.off()
eps=.10
#tiff("fig4.tiff", units="in", width=6, height=8, res=300)
par(mfrow=c(4,3));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
TriMn=RowTriMeans(IgG_raw$NSI)
ResMat=(IgG_raw$NSI)-matrix(rep(TriMn,12),ncol=12)
SubMat=abs(ResMat)<quantile(as.vector(abs(ResMat)),prob=1-eps)
for(j in 1:12){
plot(TriMn,ResMat[,j],
xlab="Raw Signal Tri-Mean",ylab="Residual Raw Signal",
main=paste("IgG Raw Signal ; Sample",j),
cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
y=ResMat[,j]
x=TriMn
fitj <- rq(y ~ x, tau = .5,subset=which(SubMat[,j])) #calculating the the quantile reggression fit, excluding the outliers
abline(coef(fitj),col=clrs[j])
}
#dev.off() ### here we are plotting the autoantigen rowmeans as a function of their residuals.
IgG_raw$NSI <- as.matrix(IgG_NSI)
IgG_raw$NSI=IgG_raw$NSI*IgG_raw$SNR
par(mfrow=c(1,1));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
plot(density(log10(IgG_raw$NSI[,1]+ 0.5)), main=" log10(NSI +0.5)", sub=" . ",type="l")
for(j in 2:12) {
dens<-density(log10(IgG_raw$NSI[,j]+0.5))
lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
IgG_raw$X_norm_1=IgG_raw$NSI
eps=.10
IgG_raw$QregFits=array(NA,dim=c(12,70))
TriMn=RowTriMeans(IgG_raw$X_norm_1) ### rowmeans
ResMat=IgG_raw$X_norm_1-matrix(rep(TriMn,12),ncol=12)
SubMat=abs(ResMat)<quantile(as.vector(abs(ResMat)),prob=1-eps)
#tiff("fig5_.tiff", units="in", width=7, height=6.5, res=300)
par(mfrow=c(2,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
plot(rep(TriMn,12),as.vector(ResMat),
xlim=quantile(TriMn,prob=c(0,.8)),
ylim=quantile(as.vector(ResMat),prob=c(0.05,.95)),
type="n",xlab="Raw Signal Tri-Mean",ylab="Raw Signal",
main="IgG Raw Signal")
for(j in 1:12){
points(TriMn,ResMat[,j],cex=0.5,pch=pchs[j],col=clrs[j])
y=ResMat[,j]
x=TriMn
fitj <- rq(y ~ x, tau = .5,subset=which(SubMat[,j]))
abline(coef(fitj),col=clrs[j])
IgG_raw$QregFits[j,]=coef(fitj)
prdct=coef(fitj)[1]+coef(fitj)[2]*TriMn
IgG_raw$X_norm_1[,j]=IgG_raw$X_norm_1[,j]-prdct #model fit
}
AAplotDens(IgG_raw,SampleIndices=1:12,clrs=clrs,useMat="norm",
xlims=NA,returnModes=T,main="Normalized IgG Spot Signal")
## [1] 4256.431 4327.828 4185.033 4327.828 3970.842 3899.444 4185.033 2614.294
## [9] 3613.856 3828.047 3114.075 3899.444
legend(2000000,0.0006,lty=1,lwd=3,col=clrs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)
plot(density((IgG_raw$QregFits[1,]),type="l"), main="(QregFits)", sub=" . ")
for(j in 2:12) {
dens<-density((IgG_raw$QregFits[j,]))
lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
plot(density(log(IgG_raw$X_norm_1[,1]+max(IgG_raw$X_norm_1[,1]))-min(log(IgG_raw$X_norm_1[,1]))), main="(X_norm_1) ", sub=" . ",type="l")
for(j in 2:12) {
dens<-density(log(IgG_raw$X_norm_1[,j]+max(IgG_raw$X_norm_1[,1]))-min(log(IgG_raw$X_norm_1[,1])))
lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
legend(400000, 0.00005,lty=1,lwd=3,col=clrs[c(1,4,7,11)],pch=pchs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)
#dev.off()
#tiff("fig6.tiff", units="in", width=7, height=6, res=300)
par(mfrow=c(2,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
eps=0.10
q.eps=0.10
lambda=.10 # incremental change
nitr=95 # number of iterations (going upto 95, 46th iteration was found to be optimum)
sad=rep(0,nitr) # sum of absolute deviations
# re-init
IgG_raw$X_1step_norm_1=IgG_raw$X_norm_1
IgG_raw$X_norm_1=IgG_raw$NSI
for(i in 1:nitr){
# snoop and grab invariants
X=IgG_raw$X_norm_1
iWRS=function(i) t.test(X[i,1:6],X[i,7:12])$p.value
qMat=matrix(rep(mapply(iWRS,1:70),12),ncol=12)
qMat[which(is.na(qMat[,1])==T),] <- 0
qcut=quantile(qMat[,1],prob=q.eps)
TriMn=RowTriMeans(IgG_raw$X_norm_1)
ResMat=IgG_raw$X_norm_1-matrix(rep(TriMn,12),ncol=12)
SubMat=abs(ResMat)<quantile(as.vector(abs(ResMat)),prob=1-eps)
SubMat[qMat<qcut]=FALSE
for(j in 1:12){
y=ResMat[,j]
x=TriMn
fitj <- rq(y ~ x, tau = .5,subset=which(SubMat[,j]))
prdct=coef(fitj)[1]+coef(fitj)[2]*TriMn
IgG_raw$X_norm_1[,j]=IgG_raw$X_norm_1[,j]-lambda*prdct #calculating actual values
sad[i]=sad[i]+sum(abs(2*lambda*prdct))
}
} #### the lowest sad or "sum of deviations" occurs at the 46th iteration.
#qcut at this i=46 is 0.0054
plot(1:nitr,sad,main="iterative adj IgG")
smoothScatter(IgG_raw$X_1step_norm_1,IgG_raw$X_norm_1, nbin = 300)
abline(0,1,col="steelblue")
plot(density((IgG_raw$X_1step_norm_1[,1]),type="l"), main="(X_1step_norm_1) ", sub=" . ")
for(j in 2:12) {
dens<-density((IgG_raw$X_1step_norm_1[,j]))
lines(dens, cex=1.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
legend(400000, 0.00005,lty=1,lwd=3,col=clrs[c(1,4,7,11)],pch=pchs[c(1,4,7,11)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)
plot(density((IgG_raw$X_norm_1[,1]),type="l"), main="(X_norm_1) ", sub=" . ")
for(j in 2:12) {
dens<-density((IgG_raw$X_norm_1[,j]))
lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
legend(4000,0.0009,lty=1,lwd=2,col=clrs[c(1,4,7,10)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)
#dev.off()
deltaG=min(as.vector(IgG_raw$X_norm_1))
IgG_raw$W_1=log(IgG_raw$X_norm_1-deltaG+1)
IgG_raw$W_2=log(IgG_raw$X_norm_1-deltaG+1)
#tiff("fig7.tiff", units="in", width=8, height=7, res=300)
par(mfrow=c(2,2));par(mgp=c(1.5,.5,0));par(mar=c(2.75,2.75,1.5,0.25))
AAModesIgG=AAplotDens(IgG_raw,SampleIndices=1:12,clrs=clrs,useMat="resid",xlims=NA,returnModes=T,main="IgG Residual (Stage 1) log Spot Signal")
legend(-0.6,14,lty=1,lwd=2,col=clrs[c(1,4,7,10)],legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1) #did not plot
# get resids..
W_1=IgG_raw$W_1
RW=RowTriMeans(W_1)
# get residual matrices
RW_MAT=matrix(rep(RW,dim(W_1)[2]),ncol=dim(W_1)[2])
W_1=W_1-RW_MAT #value - trimean
W=W_1
bw.est=unlist(apply(W,2,r.bw)) ####bivariate midvariance is something like pearsons r2. variances across columns
# now, rescale
for(j in 1:12){
W_1[,j]=W_1[,j]/sqrt(bw.est[j]) #sq. root of bivariance is an estimate of sd
}
# now, get new values..
IgG_raw$W_1=RW_MAT+W_1
IgG_raw$W_2=RW_MAT+W_1
AAplotDens(IgG_raw,SampleIndices=1:12,clrs=clrs,useMat='resid',
xlims=NA,returnModes=T,main="IgG Residual (Stage 2) log Spot Signal")
## [1] -0.003662881 0.201954756 0.300397177 0.346442825 -0.101311411
## [6] 0.350412277 -0.405371468 -0.214043860 -0.314074062 -0.192608817
## [11] -0.298196252 -0.137036483
legend(2.5,0.55,lty=1,lwd=2,col=clrs,legend=c(paste(rep("B",2),rep(c(8,28), each=1)), paste(rep("N",2),rep(c(8,24), each=1))),bty="n",ncol=1)
plot(density(log2(IgG_raw$W_1[,1]+65)-7.20), main="(logW_1) ", sub=" . ", type="l")
for(j in 2:12) {
dens<-density(log2(IgG_raw$W_1[,j]+65)-7.20)
lines(dens, cex=0.5,pch=pchs[j],col=clrs[j],cex.main=1)
}
AAplotByMean(IgG_raw,SampleIndices=1:12,clrs=clrs,pchs=pchs,MeanSubset=1:12,useMat="log",main="IgG log normalized")
AAplotByMean(IgG_raw,SampleIndices=1:12,clrs=clrs,pchs=pchs,MeanSubset=1:12,useMat="resid",main="IgG residuals log normalized")
#dev.off()